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ABSTRACT 

Context. Comets and chondrites show non-monotonic behaviour of their Deuterium to Hydrogen (D/H) ratio as a function of their for¬ 
mation location from the Sun. This is difficult to explain with a classical protoplanetary disk model that has a decreasing temperature 
structure with radius from the Sun. 

Aims. We want to understand if a protoplanetary disc with a dead zone, a region of zero or low turbulence, can explain the measured 
D/H values in comets and chondrites. 

Methods. We use time snapshots of a vertically layered disk model with turbulent surface layers and a dead zone at the midplane. 
The disc has a non-monotonic temperature structure due to increased heating from self-gravity in the outer parts of the dead zone. We 
couple this to a D/H ratio evolution model in order to quantify the effect of such thermal profiles on D/H enrichment in the nebula. 
Results. We find that the local temperature peak in the disk can explain the diversity in the D/H ratios of different chondritic families. 
This disk temperature profile leads to a non-monotonic D/H enrichment evolution, allowing these families to acquire their different 
D/H values while forming in close proximity. The formation order we infer for these families is compatible with that inferred from 
their water abundances. However, we find that even for very young disks, the thermal profile reversal is too close to the Sun to be 
relevant for comets. 

Key words, protoplanetary disks - astrochemistry - Meteorites, meteors, meteoroids - Comets: general - Planets and satellites: 
composition 


1. Introduction 

The Deuterium-to-Hydrogen (D/H) ratio is one of the most dis¬ 
cussed isotopic ratios in planetary sciences. This interest is due 
to the strong dependence it has on the temperature of an icy 
body’s formation location ( [Ceccarelli et al.|2014| >. In the inter¬ 
stellar medium (ISM), the water D/H r atio is measured to be up 
to ~ 9 x 1 CL 4 (Brown & Millar 1989), although the values are 
much higher in the organic matter of the interste llar hot cores 
(~ 1 x 10” 3 ) and cold clouds (up to ~ 2 x 10 2 ) (Dartois et al. 
2003 1 |Robert||2006[ |Parise et al.||2012| . When a protoplanetary 
disk forms around a young star, these ices are heated to more 
than few hundred K. At these high temperatures, gaseous HDO 
reacts with nebular H 2 to form water vapor and HD, 


HDO + H 2 +± H 2 0 + HD 


(1) 


and this results in a decrease of the D/H ratio in water because 
[HD]/[H 2 \ is lower than D/H in water, even accounting for 
chemical affinity (Lecluse & Robert||1994] l. The rate of this re¬ 
action depends strongly on the disk temperature. Therefore, ices 
forming in regions with higher past temperatures should have a 
lower D/H ratio. D/H measurements in Jupiter and Saturn find 
D/H ~ 2 x 10~ 5 (Lellouch et al. 2001}, the value also found in 
the ISM Hydrogen gas. This value is usually interpreted as the 
D/H ratio in the protosolar nebula’s H 2 . For the rest of this work 


we follow [Drouart et al.| ( |1999| and others in defining / as the 
D/H enrichment factor in water: 


/ = 


[■ HDO]/[H 2 0] 


( 2 ) 


[HD]/\H 2 ] 

Thus for the nebular H 2 gas, by definition, / is equal to unity. 

1.1. D/H Ratio Observations in the Solar System 

The water D/H ratio of small bodies in the solar system has been 
measured for a large number of comets fMumma & Chamley| 
201 l|l, meteorites ([Robert 2006) and in the plumes of Saturn’s 


moon Enceladus ( |Waite et aL 2009 Spence ret al.|2009) . In me¬ 
teorites, both organic materials and water contribute to the bulk 
D/H value. Therefore, separating the two components is crucial 
to understanding the contribution from water. This was done by 
[Alexander et al.| ( |2012j l, who found / values ranging from 3 to 
5 in carbonaceous Cl chondrites, and ~ 15 in ordinary chon¬ 
drites (OCs). More recently, a higher value reaching f ~ 100 has 
been also inferred for OCs ( jPiani et al.|2015| ). This implies that 
either low D/H minerals still coexist in the matrix mixed with 
deuterium-rich organics, or an unknown process has fractionated 
water and organics differently from an initial low D/H reservoir, 
before their incorporation in the matrix. In neither case the low 
values for the D/H ratio obtained from this method have ever ex- 
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Table 1: The measured water D/H fractionation and abundance in small bodies in the solar system. 


Source 

/ 

% Water Abundance 

Ordinary Chondrites (OCs) 

LL3 

29-44 

- 

Semarkona (LL3) 

-100 

- 

H and L 

< 15 

1 

Carbonaceous Chondrites 

Cl 

3-5 

20 

CM 

-4.5 

13 

CR 

- 8 

6 

CO 

4-7 

1 

CV 

1-4 

1 

Comets 

103P/Hartley 

-7.4 

- 

67P 

-21 


Oort Cloud Comets (OCCs) 

10-25 

- 

Moons 

Enceladus 

- 15 

40-50 


isted in chondrites. Taken at face value, these results contradict 
the predictions from the expected formation locations, since CIs 
are usually associated with C-complex asteroids that formed be¬ 
tween the giant planets, and OCs are associated with S-complex 


JFC 67P/Churyumov-Gerasimenko by Rosetta, where the value 
found lies well within the OCCs range (Altwegg et al.|20l5 i. A 
summary of all of the D/H ratio observations discussed here, and 
the water abundances in chondrites, can be found in Table 1. 


asteroids that formed in situ sunward of Jupiter’s orbit (Chap 
man|1996 Mothe-Diniz et al.|2003[ [Walsh et al.|201 l) . |Alexan- 


der et al. ( 201 0| > noted though that due to Hydrogen escape via D/H Ratio Connection to Protoplanetary Disc Models 


Fe oxidation, the D/H ratio in OCs should be treated as an upper 
value and could have originally been much lower. This effect was 
found to be limited to OCs, allowing us to use the values mea¬ 
sured in other chondritic families to understand the processes at 
play. 

The carbonaceous chondrites, CRs, are thought to have a 
high D/H value ([Robert & Epstein|[T982[ l. The value measured 
by [Alexander et alT ( |2012| ) of / ~ 8 is very surprising, since this 
is almost twice the value in CIs, and is actually higher than that 
found in the comet 103P/Hartley. The D/H ratio was measured 
also in the ordinary chondrites type LL3 and found to be even 
higher than that in L and H type ordinary chondrites ( |Robert| 
|2006[ [Alexander et al.|2010| l. Measurements of the D/H ratio in 
the carbonaceous chondrite types CO, CM and CV were found 
to be also low providing further evidence for this discrepancy 
in the carbonaceous and ordinary chondrite D/H ratios (Robert 
|2006[[Alexander et al.|2012] ). 

In Oort cloud comets (OCCs), classically thought to have 
formed between the giant planets, / was found to range between 
10 and 25. On the other hand, the D/H ratio was measured for 
the first time in the Jupiter family comet (JFC) 103P/Hartley, 
a family classically thought to have formed in the Kuiper belt 
beyond Neptune, and was found to be ~ 7.4, the value also 
found in Earth’s oceans (Hartogh et al. 2011 1 . This was sur¬ 
prising since JFCs are thought to have formed in an area fur¬ 
ther out than the OCCs (|Morbidelir| |2008j l (although cf. the 
model of [Brasser & Morbidellil|2013[ where all comets form 
beyond the orbital radius of Neptune), their D/H ratio was pre¬ 


dicted to be higher than or at least equal to the OCCs (Mousis 
|et al.||2000[ |Kavelaars et al.||2011] >. What made the situation 
even murkier is the recent measurement of the D/H ratio in the 
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Classical disk models with a standard monotonically decreasing 
temperature profile ( |Lynden-Bell & Pringle|1974[ Pringle|1981| 
all predict that the D/H ratio should similarly follow a mono¬ 
tonic profile, increasing with the icy body’s formation distance 
from the sun ( |Drouart et al.||1999[ |Mousis et al.||2000] |Hersant| 


et al.||200T] l. An alternative model is that of Yang et al. (2013 1 
who use a new D/H 2D model including the infalling material 
from the cloud, giving a non monotonic D/H ratio profile in the 
nebula due to constant influx of unequilibrated water. However, 
they did not discuss the D/H ratio in chondrites, and their mono¬ 
tonic temperature profile out to a radius of 10 AU cannot explain 
the diversity of D/H ratios found in the inner solar system. An¬ 
other model is that of [Jacquet & Robert] ( |2013| , who tried to 
explain the chondritic diversity with a classical disc model that 
includes an interplay of inward advection and outward diffusion 
in the nebula. This model however also predicted a monotonic 
D/H profile, and although it can broadly explain the chondrites 
D/H range, it does not explain for example why CRs, that formed 
closer to the sun than CIs ( jWood|2005| ), have a higher D/H ratio. 

The observed D/H fractionation variations found in chon¬ 
drites are a challenge to classical disk models, since the parent 
bodies of most chondrites should have formed in the same gen¬ 
eral region, except those of CIs that probably formed few AU 
further out in the disk. All these results indicate that either the 
D/H evolution models used are incomplete or that the thermal 
profile in the protosolar nebula was not monotonic, the hypothe¬ 
sis we are going to explore in this work. 

Turbulence within protoplanetary disks drives outward angu¬ 
lar momentum transport that allows material to spiral in and be 
accreted onto the forming star (e.g. Pringle|[T981 ] >. This turbu¬ 
lence is thought to be driven by the magneto-rotational instabil- 
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ity (MRI) (Balbus & Hawley|l99l) l. However, it is now generally 
accepted that the midplanes of protoplanetary disks have a region 


We discuss the uncertainties associated with the a parameter in 
Section 3.5. 


of zero or weak turbulence known as a “dead zone” (e.g. Gam 


mie 1996; Fromang et al. 2002 

Martin et al. 2012a, Dzyurkevich 

et a .|2013} Turner et al. 2014; 

Cleeves et al. 2015 ). This leads to 


We follow Martin & Lubow (2011) and model the protoplane¬ 
tary disc as a layered accretion disc (see their equations 1-16). 
The surface density evolves according to conservation of mass 
and angular momentum (jPringle||~1981 1 . The temperature struc¬ 


ture is governed by a simplified energy equation that balances 
viscous heating with black body cooling (see e.g. Pringle et al. 


(1986|>; Canizzo (1993])). There are two types of turbulence in 


protoplanetary discs, magnetic turbulence driven by the MRI and 
gravitational turbulence ( Paczynski||1978| ). 

The MRI requires a critical level of ionisation to operate so 
that the gas is well coupled to the magnetic field. This can be 
achieved if the temperature of the disc is larger than the critical 
T > T cr i t = 800 K ( |Umebayashi|1983~] ). In this case, the MRI op¬ 
erates at all disc heights. However, the temperature in the outer 
parts of the disc is lower than this. In this region, the disc surface 
layers may be ionised by external sources of ionisation such as 
cosmic rays or X-rays from the central star to a_maximum sur¬ 
face density depth of X c 


200gem 2 (e.g. 


Gammie 


(1996); 


|Glassgold, Najita & Igea| ( |2004[ )). If the total surface density of 
the disc is larger than this, X > X crlt , then there is a dead zone at 
the midplane with surface density Xd = X - X cri | and the active 
layers have X a = X cr j t . Otherwise, if X < X cnt , then there is no 
dead zone layer so that Xd = 0 and X a = Z. 

The MRI active layers have a Sha kura & Sunyaev] (|1973 i 
viscosity a parameter of 0.01 (e.g. Hartmann et al.[( |l998j )). In 


a gradual accumulation of gas in the dead zone region, resulting 
in an increase in the temperature and pressure (jArmitage et al.l 
[2001] |Zhu et al.|[20T0l |Martin & Lubow||2011| |2013) . The in- 
creased surface density of the disc can lead to a second type of 
turbulence, driven by gravitational instability ( ]Paczynski]|1978[ 
|Lodato & Rice|2004l ) and this results in an increase in the tem¬ 
perature locally in that region. The question we are going to 
tackle in this work is: What effect does a local peak in the tem¬ 
perature have on the D/H ratio in protoplanetary disks and can 
it resolve the discrepancies in the observations of the D/H ratio 
in small bodies in our solar system? In section 2 we discuss the 
model we use to simulate the processes involved. In section 3 we 
discuss the results and implications. We draw our conclusions in 
section 4. 


2. The numerical model 

2.1. The protoplanetary disk model 


Gravitational turbulence requires the Toomre ( 1964| parame¬ 
ter to be less than the critical, Q < gerft = 2. While a dead zone is 
present in a protoplanetary disc, the flow through the disc is not 
steady because material accumulates in the dead zone. With suf¬ 
ficient material in the dead zone, it may become self gravitating, 
thus a small amount of turbulence may be driven. We include 
an additional viscous term in the surface density evolution equa¬ 
tion and a heating term in the temperature equation. This extra 
heating in the massive dead zone can eventually cause the disc 
to reach the critical temperature required for the MRI. Once this 
is reached, there is a snow plough effect through the disc and the 


whole disc becomes MRI active in an outburst phase (Martin & 
|Lubow|2013] ). As material drains on the star, the disc cools and 
the dead zone can reform causing repeating outburst and quies¬ 
cent cycles. Once the infall accretion rate onto the disc drops off, 
there may not be sufficient inflow through the disc for another 
outburst to occur, but there can still be a dead zone within the 
disc. 

The infall accretion rate onto a forming star varies in time 
although the details of the evolution depends on the specific disc 
model (e.g. Shu||1977| Basu|1998] |Bate|201 1| |Kratter, Matzner 


our model, the dead zone has zero turbulence, unless it becomes 
self-gravitating. Turbulence in the dead zone could be driven by 
other sources such as hydrodynamic instabilities (including the 
baroclininc insta bility, |Klahr & Bo denheim er|(2003| ); [Petersen, | 
Julien & Stewart (20071; Lesur & Ogilvie ( 2010| i) or be induced 


from the magnetohydrodynamic instability in the active surface 
layers ({Fleming & Stone]|2003| |Simon, Armitage & Beckwith] 
|201 1 [[Gressel, Nelson & Turner|2012| ). However, we note that a 

small amount of turbulence within the dead zone does not signif¬ 

icantly alter the qualitative disk structure and behaviour (Martin 

& Lubow||2014)l. The conclusions of our work are not signifi- 


j& Krumholz|20()8| 7However, it is thought that at early times the 

infall accretion rate is approximately A7 m[a ii = c\/G. Thus, as¬ 
suming a cloud temperature T ~ 10K, the initial infall accretion 
rate is ~ 10 * 5 M G yr 1 . In this work, we consider several con¬ 
stant infall rates and thus analyze the disc structure at different 
evolutionary times. We choose three infall accretion rates onto 
the disc of Mj n f a u = 2 x 1CT 5 , 1 x 10~ 6 and 1 x lO~ 8 M 0 yr _I . 

We solve the accretion disc equations with 200 grid points 
spaced equally in log R. For each infall accretion rate, we run the 
model until the disc reaches either a steady outburst cycle (for the 
higher two accretion rates), or a steady state (fully MRI active) 
disc solution (for the lowest accretion rate). These models repre¬ 
sent different stages of the disc evolution. In Fig. |T] we show the 
surface density and temperature profiles for the disc as a function 
of radius at three different times for M m f a n = 2 x 10 5 M G: yr 1 . 
The first is a time immediately before an outburst, the second 
is during the outburst, and the third is immediately after the out¬ 
burst. During this early phase of the disc evolution, the timescale 
between the outbursts is only 3 X 10 3 yr. As we expect the accre¬ 
tion rate on to the disc to decrease exponentially in time, the 
timescale between the outbursts increases during the disc life¬ 
time (see e.g. Martin et al.| |2012b) . Furthermore, the radius of 
the temperature peak (the radius at which the disc becomes self- 
gravitating) moves inwards in time (Martin & Lubow [2013) ). In 
Fig. [2] we show the surface density and temperature profiles cor¬ 
responding to the different disk stages and infall accretion rates. 
For high accretion rates, the local thermal peak is present in the 
disk, and disappears for M m ia\\ ~ 1 x 10 8 M 0 yr 1 because there 
is no dead zone at this low accretion rate. The model behaves 
as a classical accretion disk without a dead zone (see for ex¬ 
ample for comparison jBaillie et al.||2015~| who found tempera¬ 
ture and surface density profiles comparable to the model we 
use for a similar accretion rate value) . The position of the ther¬ 
mal peak is also sensitive to the infall accretion rate. The model 
with Mi n f a ii = 1 x 10 6 M G yr 1 is used as the nominal disk model 


candy affected by an additional source of turbulence within the 
dead zone, unless the source is strong enough to be able to pro¬ 
duce a steady state disk. This would require a level of turbulence 
comparable to that produced by MRI (Martin & Lubow||2014|>. 


in this work, and this corresponds to the model used in (Martin 
|& Lubow] 2013} ). We note that there is some uncertainty in the 
value of Z cr it such that it could be much smaller, which would 
mean that the temperature peak would exist for even lower ac¬ 
cretion rates (see section 3.5 for more on this point). 

Article number, page 3 of|9] 



















































































































A&A proofs: manuscript no. dtohdeadzone 




Fig. 1: The disk surface density (left panel) and temperature (right panel) structure for the high infall accretion rate M; n f a u = 
2 x 10 5 M 3 yr 1 . The profiles are shown just before, during, and just after an accretion outburst. The “during outburst” profile is 
used in the calculations corresponding to Fig. [4] 
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Fig. 2: The disk surface density (left panel) and temperature (right panel) for Mj n f a n = 2 x 10 5 M 0 yr 1 (at a time just before an 
outburst), 10 6 M 0 yr 1 (at a time between outbursts), and Mi n f a ii = 10 s ;V/ 3 yr 1 (in steady state). The thermal gradient reversal 
observable in the right panel is main element of the disk model. We use M; n f a ii = 1 x 10 6 M e yr -1 as our nominal disk model. 


2.2. The D/H evolution model 


diffusivity is 


Dust and grains in the protoplanetary disk settle to the disc mi- 
plane within the dead zone layer as there is little or no turbulence 
there. Thus we expect that the solar systems bodies formed in the 
dead zone layer. We couple the protoplanetary disk profiles to a 
classical 1-D D/H ratio evolution model. The D/H ratio in the 
dead zone layer is assumed to follow 


Of 

dt 


— kP (A - f) + 


i a 

l^RdR 



df_\ (2 k dY. d 
dRl + \E d dR 



(Drouart et al.T1999), where k(T) is the rate of isotopic exchange, 
P is the gas pressure, A(T) is the fractionation at equilibrium and 
Ed is the gas surface density in the dead zone layer. The turbulent 


where Vd is gas viscosity in the dead zone layer and P, the 
Prandtl number. Vr is the radial gas velocity: 

v ‘ = -z^Tr^‘ r,I1) <5) 

(IPringle| 1981) . This equation is valid in the midplane deadzone 
layer assuming that Ej » E„, which is satisfied in the disk pro¬ 
files that we use (cf. section 3.5 for caveats discussion). The vis¬ 
cosity due to the self-gravity in the dead zone layer is given by: 
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(although we note that any value of order unity will lead to qual¬ 
itatively similar results). 

An explicit Forward-time Centered-space (FTCS) scheme 
was first used to solve this equation. During tests in the case 
of a regular monotonic thermal profile, the code converged if 
a sufficiently small time step was used. When the new ther¬ 
mal profiles (with temperature gradient reversals) were used the 
code diverged even for very small time steps, except for almost 
vanishing turbulence. Thus, we employ semi-implicit (Crank- 
Nicholson) and fully implicit schemes to solve the equation over 
the same grid used in the disk model, but we take the water 
snow line (the radial location in the disc inside of which water 
is gaseous, and outside it is solid, that occurs at a temperature of 
around 7' snow = 170 K |Lecar|2006| > as the outer boundary, since 
the deuterium exchange can occur only in the vapor phase. 


Fig. 3: The dead zone midplane turbulent viscosity parameter a^. 
If the disc is not self-graviting (i.e. Q > Q CI \ t ), then = 0. Oth¬ 
erwise, its value is depicted by equation (|7j. Note that we manu¬ 
ally smooth the viscosity variation (hence ad) over 5 grid points 
to avoid code divergence. For M; n f a n = 10 6 M 0 yr 1 the smooth¬ 
ing is between 2.2 and 2.5 AU and for M m Mi — 2 x 10~ 5 M Q yr -1 
this is between 5.2 and 5.5 AU. 


C * 2 * * 

Vd=a d -^-, (6) 

where a,j is the turbulence parameter. For all this work, a a, X,/ 
and v,i are the parameters in the midplane (MRI inactive dead 
zone), and not the value in the MRI active disk surface. C s is 
the speed of sound and O the gas Keplerian velocity. The sound 
speed is calculated with the midplane disc temperature profile. 
In this disk model, a d is zero unless the disc locally satisfies 
Q < 2crit- and then we have 


a d 



if Q < Gcrit 
otherwise 


(7) 


([Martin & Lubow||201 lj ). We note that if we were to include 
an additional source of visosity within the dead zone, other than 
self-gravity, then we would add an additional term to this equa¬ 
tion (see |Martin & Lubow|2014[ ). Fig.[3]shows the viscosity pa¬ 
rameter in the dead zone for the highest two disk infall accretion 
rates that we consider. This viscosity in the dead zone is much 
smaller than that in the active layers. 

The first term on the right side of eq. 3 describes the chemical 
isotopic exchange between HDO and H 2 , and takes the pressure 
and temperature of the disk as the input. The second term de¬ 
scribes the gas diffusion due to turbulence that depends on a and 
X. The third term describes the gas advection that is dependent 
on a, X and Vr. In contrast with previous models, we do not ne¬ 
glect this term and we solve the entire equation. The effect of a 
temperature peak in the nebula in our model is thus controlled 
by the interplay of these terms. In the absence of turbulence, the 
first term dominates and the effect of the temperature peak is to 
decrease the local D/H ratio. However, if the turbulence is strong 
enough, the diffusion and advection terms erase this gradient. We 
run all simulations with k(T) and A(T) from (Lecluse & Robert 
|1994[ >, and set P r to 0.5 ( |Prinn|[T990l |DubruI e & Frisch| 1991 1 > 


3. Results & discussions 


Most simulations in the literature begin with a spatially constant 
/ = 25, which is close to the highest value observed today in 
the solar system of 29—44 in LL3 meteorites (except Semarkona 
where a possible value of up to 100 was recently inferred 
|ani et al.1[2015| )). However, since we are interested in the D/H 

ratio difference between two bodies rather than the absolute val¬ 
ues, we begin our standard simulations with /o = 15 which is 
the average value found in comets. The model evolves from this 
value and we check the effect of the thermal gradient reversal 
on the D/H ratio profile. We are therefore implicitly assuming 
that another transient heating process decreases the initially very 
high (LL3 or even interstellar) D/H ratios to the lower values we 
are using as the initial condition. A possible process for this is 
the gravo-magneto disc instability (Martin & Lubow pOl 1 1 and 
its associated accretion outburst that we first test here (see also 
Owen & Jacque fp015) >. The accretion outburst occurs when the 
local peak in the temperature profile becomes high enough that it 
reaches the critical temperature, T cut , required to trigger the MRI 
in the dead zone. During the outburst the disc becomes MRI ac¬ 
tive throughout and a large amount of material is accreted onto 
the Sun in a short time. After the outburst, the disk cools, the 
dead zone reforms and providing that there is sufficient accre¬ 
tion inflow, the cycle repeats. 


We first run a test simulation with a disk profile representa¬ 
tive of the conditions during an accretion outburst (for M; n f a u ~ 

2 x 10 5 M Q yr cf. Fig. m). We begin with an enrichment fac¬ 

tor of f) = 35 (representing the very high D/H enrichment that 

can be found in LL3 or certain ISM environments). The viscos¬ 

ity parameter is a — 10 2 everywhere. Further, this leads to fast 
transport of material. Thus, there is very fast D/H ratio evolu¬ 
tion as shown in Fig. [4] where cometary values are reached in 
the inner 10 AU over the timescale of an outburst. The outburst 
period lasts a few hundred to few thousand years and during this 
time the D/H enrichment decreases to cometary or lower values 
in the inner disk, but remains high in the outer disk. These out¬ 
bursts can happen multiple times during the disk phase and they 
alter the D/H ratio inhomogeneities that may have existed prior 
to the instability. For any process to have measurable effect on 
the D/H ratio today it has to happen after the last accretion out¬ 
burst, when the disk infall rate from the cloud has dropped to 
values below certain threshold point. Thus, in the next Section 
and for our standard model we consider the evolution at lower 
infall accretion rates. 
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3.1. Standard Model 

We now describe our standard model with an infall accretion rate 
of Minfaii = 10 6 M 0 yr 1 . We consider a time between outbursts 
as shown in Fig. [2] We assume that the last accretion outburst 
happened at an infall accretion rate of around Minfaii reached 1 x 
10 6 M 0 yr 1 . Thus, the disk structure does not change rapidly 
after this. The dead zone accumulates material and heats up in 
the outer parts by self-gravity leading to the thermal gradient 
reversal, but without reaching sufficiently high temperature (~ 
800 K) to trigger a further outburst. Results for this model are 
shown in Fig. [5] We notice that the D/H ratio is decreasing in 
mainly two locations: the inner hot disk, and around the thermal 
peak centered at 3.5 AU. The width and limits of the D/H ratio 
dip around this region are controlled by the turbulence strength. 
We notice also that in the inner disk, equilibrium is reached in 
less than 10 4 yr, faster than the evolution time of the disk, thus 
justifying the use of a snapshot disk model in this region. We 
discuss this further in Section 3.5. Finally, it should be noted 
that Fig. [5] shows / decreasing to values lower than the average 
cometary D/H in the outer disk at late times, but this is an artifact 
of our disk cooling handling as discussed in the next section. 



Fig. 4: The D/H ratio evolution during an accretion outburst. This 
simulation starts with / = 35 as initial condition to reflect the 
very high enrichment values found often in LL3 and certain ISM 
regions. 


3.2. Disk cooling and photoevaporation 

As long as there is a hot region (T higher than ~ 400 K) along 
with turbulence within the disk, the global D/H ratio will con¬ 
tinue to evolve until it reaches A(T) (/ ~ 1) throughout the disk 
(all water equilibrated to nebular gas D/H ratio). Such values of 
/ for the entire disk are contradictory to observations. There are 
two possible explanations for this. 

First, the disk could have cooled down sufficiently quickly 
that the D/H ratio profile became frozen, as considered in the 
classical models ( [Mousis et al. ||2000} Hers ant et al.i 2001) . For 
low enough temperatures, the chemical exchange (through k(T)) 
becomes very slow and inefficient. The gas-gas reaction then 
stops completely once water condensed into ice. This needs to 
happen before the value of / in the chondrites region becomes 
too low. This is equivalent to the water snow line radius mov¬ 
ing quickly inside the chondrites formation region. The snow 
line is the radial location in the disk inside of which water is 
gaseous, and outside it is solid, that occurs at a temperature of 
around 7/ now = 170 K ( |Lecar|[2006j ). For our standard model 
with an accretion rate of M m f a n = 1 x 10 -6 M 0 yr -1 , the snow 
line is at a radius of around 9 AU. The disk will remain in this 
state while M; n f a n decreases in time and the snow line moves 
inwards slowly ( [Martin & Livio| ( [2012| ) and cf. equation 19 in 
Martin & Livio;( 2013| >). For the disk parameters we have cho¬ 
sen, our model shows that when M; n f a n reaches ~ 10 s M 0 yr -1 , 
the disk will quickly become cool with a classical monotonic 
thermal profile with a snow line radius of around 1 AU (see Fig. 
[2] bottom panel). The infall accretion rate is given by: 


3/n tall 


= M,exp(-l) 


( 8 ) 


(equation 19 in 


Martin et al. 


2012b), where M\ is the initial in¬ 
fall accretion rate (in this case we take M\ = 10 -6 M 0 yr -1 ), t is 
time, f ff is the free fall time scale (we take f(f = 10 5 yr e.g. 


10 - 


Ar- 


1 AU, thus freezing the D/H ratio in the chondrites regions. The 
final D/H profile in Fig. [5]is thus in the solid phase. 

The second explanation for how the D/H radio profile be¬ 
came frozen is that the disk was completely photoevaporated on 
a similar timescale to the evolution of/. During the photoevapor- 
tation process at the end of the disk li fetime the disk is dispersed 


mitage et al. 20011). This equation shows that Mi n f a ii will reach 
8 M 0 yr -1 in about 5x10 s yr. We will hence stop our sim¬ 


ulations at this time, and use the end state results to fit the mea¬ 
surements. We are therefore implicitly assuming that the transi¬ 
tion to a cold classical MRI active disk happens quickly, with 
the snowline moving in from the giant planets region to around 


on a short timescale of around 10 s yr ( [Clarke et al.|2001[[Alexan¬ 
der et al.|2006| Owen et al.|2010 i. Thus, the profile for the ratio 
of D/H would become fixed in the chrondrites at this time and 
similarly, the shape of the profile shown in Fig.[2]would become 
fixed at this time. 


3.3. Implications for chondrites 

Now we discuss the implications of our model for chondrites. 
Figure [5] shows our nominal D/H profile obtained for our nom¬ 
inal model, with some known chondritic D/H ratio ranges. The 
profile shows a clear peak reaching / ~ 15 around 2 AU, fol¬ 
lowed by a dip reaching / ~ 4 around 3.5 AU. The curve then 
increases all the way to the cometary values. The absolute peak 
and dip positions with respect to the sun in this plot are not too 
important, since these are sensitive to and thus can be var¬ 

ied (since the chondrites parent bodies probably formed slightly 
outward of these positions, to allow for S-type asteroids to form 
around 2.2 AU). We chose a single model for simplicity, and 
this particular value to remain consistent with Martin & Lubow 
( 2013] ) to allow for comparison. As seen in Figure [5] this pro¬ 
file can fit the D/H values found in most chondrites, and explain 
their diversity from their formation location. Thus, this explains 
why CRs have a higher D/H value than other type of chondrites. 
The relative formation locations in this plot are also compati¬ 
ble with their water abundances ( jWood|[2005[ |Brearley [2006), 
although this abundance alone has its limitations due to the very 
nature of chondrites as an association of different components. 
Other indicators such as chondrules abundance should be taken 
into account. CIs and CMs are more aqueously altered than CRs 
and OCs that are more reduced. Assuming naively that water 
abundance increases with heliocentric distance (although cf. the 
dynamical water distribution models of for example |Ciesla &' 


|Cuzzi| ( |2006| ); |Ali-Dib et al.| ( |2014| l), this implies that CIs and 
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Fig. 5: The D/H ratio enrichment evolution for our nominal disk. 
These profiles allow an important D/H ratio diversity in the chon¬ 
drites formation region. They explain the values measured in 
the different chondrites families from their formation location. 
These inferred formation locations are also compatible with the 
water abundances of chondrites, with OCs and CRs containing 
less water (hence probably forming closer to the sun) than CMs 
and CIs. The exact distance relative to the sun is not important 
since the model is scalable with Minfaii- The D/H ratios in comets, 
Enceladus, and VSMOW (Earth) are shown for reference. Hav¬ 
ing / decreasing to low values in the cometary region at late 
times is an artifact of our disk cooling handling. It should be 
higher with consistent cooling treatment. Indeed, when tempera¬ 
ture drops below water condensation point, we keep evolving /, 
letting it continue to drop, while it should be fixed at the value 
reached when the water became isolated from gas by freezing. 


CMs formed further out than CRs and OCs, which is compatible 
with this profile, and as proposed by Wood ( 2005| >. We can ten¬ 
tatively try to fit the remaining chondrite families as well (COs 
and CVs) that are more reduced than CMs but less than OCs. The 
D/H profile does contain regions that fit their values, however we 
note that there are many uncertainties present in the model since 
it is only a proof of concept. A caveat in this model is that CIs 
form too close to CMs, although a distance of several AU might 
be needed between the formation locations of the two to allow 
the formation of chondrules. 


3.4. Implications for comets 

The next step is to check if such model can explain the D/H ra¬ 
tio diversity in comets. Naively one can expect that earlier in the 
disk lifetime, when M m Mi was higher, the thermal peak could 
have existed further out in the disk, maybe in the comets region. 
This can lead to a D/H profile analogous to Fig. [5] but further 
out in disk, with its own peak and dip. This dip can give a neat 
explanation for the relatively low D/H value in 103P/Hartley, 
with the other comets forming in other locations. Assuming a 
classical comets formation model with the JFCs forming fur¬ 
ther out than OCCs, 103P/Hartley could have formed on the 
D/H dip, while 67P/C-G formed slightly further out outside of 
the dip, and the OCCs forming inside of it. To test this hypoth¬ 


esis we used a profile derived from our disk simulation with 
Minfau = 2 x 1(T 5 M g yr 1 corresponding to the time just before 
an outburst. Results are shown in Fig. [2](top panel). At this stage 
of the disk evolution, the thermal reversal (and the correspond¬ 
ing D/H enrichment dip) are around 8 AU, further out than in our 
canonical case, and as expected from an early disk. The thermal 
reversal’s position, even for such young disk, is still too close 
in to be relevant for comets formation. Another problem posed 
by this profile is that it almost certainly leads to an outburst, ho¬ 
mogenizing the D/H ratio across the disk. Within the model and 
parameters we use, we are unable to explain the low D/H ratio in 
67P/C-G. The possibility that there exists another set of param¬ 
eters and/or assumptions within the same framework that can 
lead to a thermal reversal in the comets region is not excluded 
though, and is left to future work. We note that all cometary D/H 
models assume that the value measured in the comet’s ejecta re¬ 
flect its bulk value, although |Podolak et al.| ( |2002) l showed that 
the nuclei D/H ratio might be different that on the surface. Ad¬ 
ditionally, experiments by Brown et al.] ( j2012) > showed that the 
measured value might also change a function of the instrument- 
target distance. 

Another seemingly unrelated problem that can be addressed 
using such models is the origin of crystal silicates and CAIs in 
comets. These minerals can form only at temperatures in the 


(jCampins & Ryan 

1989 

Wooden et al. 2005; Chi et al. 

20091 

Kelley & Wooden 

2009 

). How the high temperature minerals 


got to the cold region where comets form is a classic prob¬ 
lem. Some of the proposed solutions were outward turbulent 
diffusion of particles (Bockelee-Morvan et al. 2Q02j > and pho¬ 
tophoresis (Mousis et al.; |2007j l. The recent observation of nar¬ 
row crystal silicates features in the spectra of a young solar 
like star during an accretion outburst indicated that these out¬ 
bursts might be the formation mechanism of high temperature 
minerals (Abraham et al. 2009). For the accretion outburst for 
Afinfaii ~ 2 x 10~ 5 M 0 yr , the outburst trigger radius (the radius 
of the temperature peak) in our model is around 7 AU (cf. Fig. 
[lj, considerably widening high temperature minerals formation 
region. It should be mentioned that most materials inside of the 
trigger radius get accreted onto the sun during the accretion out¬ 
burst, so only elements forming beyond this radius remain in the 
disk. Quantifying any of these possible solutions is beyond the 
scope of this work. 


3.5. Caveats 

Since this work was only intended to be a proof of principle 
highlighting the concept and quantifying the relative strengths 
of diffusion and chemistry in a non monotonic nebula, numerous 
assumptions and simplifications were made in this model: 

- Ideally, one should use a time evolving disk model coupled 
dynamically with the D/H module, to track the simultaneous 
evolution of both components. However, for simplicity we 
used a static (snapshot) disk profile with the time dependent 
D/H module. Hence we are making the implicit assumption 
that the D/H ratio evolves on a shorter timescale than the 
disk. This assumption is justified by the short timescale of the 
D/H evolution (10 4 - 10 5 yr) compared to the disk evolution 
timescale (~ 10 6 yr). 

- Our disk profiles are derived from a layered (active and dead 
zones) disk. In this work, we are only tracking the D/H ra- 
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tio in the midplane (dead zone) and ignoring any effect the 
active layer might have, including the sedimentation of equi¬ 
librated water. Our work is thus valid only if the dead zone 
surface density is much higher than the active layer surface 
density. Since our domain starts beyond the dead zone inner 
boundary at 0.5 AU and extends out to the snow line radius 
at 9 AU (for the particular choice of the infall rate and thus 
disk age), much closer than the dead zone outer boundary at 
23 AU. This validity condition for our model is thus applica¬ 
ble throughout the entire domain. 

We note that there are several unknown parameters in the 
layered disc model. For example, the critical surface density 
that is ionised by external sources is not well determined. 
Dead zone models that include more physics generally find 
active layer surface densities that may be very small (e.g. Bai 
( (20111 >). However, such small active layers cannot explain 


accretion rates observed in T Tauri stars (e.g. Perez-Becker 
|& Chiang] ( |201 1[ >; |Martin et al.| ( |2012b| l). Thus we fold all 
of the uncertainty into the parameter Z crit . The value of a in 
the active layers is also not well determined. However, these 
parameters do not affect the qualitative behavior of the disc. 

- In this work we started our main simulation (Fig. [5j from 
a constant D/H value throughout the disk. Realistically how¬ 
ever the preceding outburst may lead to a heterogeneous D/H 
distribution. Quantifying this effect needs a fully time depen¬ 
dent coupled disk-D/H evolution, for a smooth temperature 
variation to occur. This is left for future work. 

- Recently, 2D (r-z) steady-state models of protoplanetary 
discs have been constructed with an or-variation over the disc 
height to mimic the effects of a reduced (but non-zero) a in a 
dead zone ( [Bitsch et al.|2014j l. These disc models do not find 
the increase in temperature in the dead zone region present in 
our disc models because self-gravity does not operate. There 
is not a sufficient build up of material in their disc models to 
cause the disc to become self-gravitating. This is because the 
chosen values for a are high enough that a steady state disc is 
found. [Martin & Lubow| ( 2014| l showed that even with some 
turbulence in the dead zone, the qualitative disc behavior is 
as we have described in this paper, unless the turbulence in 
the dead zone is comparable to that in the active layer, where 
a steady state may be found. Previous 2D (r-z) simulations 
that are time-dependent and included a dead zone with a 
smaller viscosity agree with the numerical models used in 
this work ( |Zhu et al.[2009| ). Further detailed magnetohydro¬ 
dynamic time-dependent numerical simulations are required 
in order to determine the correct value of a in the dead zone 
(see for example Simon et al. |2013| ). 

- In our simplified model, we set the viscosity in the dead zone 
to be zero, except where it is generated by self-gravity. It is 
possible for other hydrodynamical instabilities to operate in 
the dead zone (for example the baroclinic instability, vertical 
shear instability and others ( Turner et al.|2014| . However, 
as discussed by Martin & Lubow[( j2014| l, the qualitative disc 
behavior is the same even if there is a small amount of tur¬ 
bulence in the dead zone. The temperature peak is still there, 
and the outbursts still occur. 


4. Conclusions 

We have coupled a D/H enrichment code including diffusion, ad- 
vection and chemical exchange to snapshots from protoplanetary 
disk model that includes a dead zone. The disk model contains 
a local temperature peak at a radius of around 3 AU due to the 
heating by self-gravity in the outer parts of the dead zone. We 
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found that this leads to a dip in the D/H profile around the same 
region, in contrast with the classical monotonic D/H models. The 
new profile can explain the origin of the D/H ratio variations be¬ 
tween the different chondrites families. We propose that Cl chon¬ 
drites (that have a relatively low D/H ratio) formed in the region 
of the thermal gradient reversal, but CRs (that have a high D/H 
ratio) formed just inside of this region. The new D/H profile also 
accommodates the formation of COs, CVs, and CMs. However, 
even with a younger disk profile the model is unable to explain 
the the D/H ratio in 67P/C-G. The thermal gradient reversal is 
too close to the Sun to be relevant. Finally we proposed that the 
accretion outbursts associated to these models can explain the 
presence of high temperature minerals across the disk. 

This work shows that detailed temperature profiles from 
time-dependent layered disk models provide a potential explana¬ 
tion for the rich variation of D/H ratios found in the solar system. 
A more detailed understanding of the role the thermal inversions 
in dead zones and outbursts plays in shaping the chemistry of the 
nebula necessitate a more elaborate exploration of the parame¬ 
ters space, to be the subject of future works. 
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